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1  INTRODUCTION 

In  recent  years  there  have  been  a  number  of  papers  dealing  with  methods 
for  remote  temperature  sensing  based  on  applying  transform  theory  to  the 
radiative  tranfer  equation.  King,  Hohlfeld,  and  Kilian  (1989)  have  shown 
that  using  differential  inversion  it  is  possible  to  successfully  determine  tem¬ 
perature  profiles  from  measurements  of  the  upwelling  intensities  in  the  at¬ 
mosphere.  In  this  research  project  an  alternative  transform  method  based 
on  optical  measure  theory  is  developed.  Previous  work  on  this  method  has 
been  carried  out  by  by  King  and  Leon  (1989,1990)  and  by  Leon  (1990).  The 
optical  measure  theory  method  has  the  advantage  that  it  does  not  require 
any  computations  of  numerical  derivatives.  Instead  the  radiance  profile  is 
approximated  by  a  rational  function.  The  Planck  intensity  is  then  deter¬ 
mined  as  the  inverse  transform  of  the  rational  function. 

In  Section  2  we  present  techniques  for  approximating  the  radiance  pro¬ 
file  based  op  rational  interpolation  and  nonlinear  least  squares  data  fitting. 
Optical  measure  theory  permits  the  inversion  of  non-Laplacian  exponential 
kernels.  Consequently,  in  Section  3,  we  extend  the  mathematical  models  to 
include  generalized  exponential  weighting  functions.  In  Section  4  we  discuss 
alternative  nonlinear  data  fitting  techniques.  In  particular  we  introduce  the 
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concept  of  total  nonlinear  least  squares  data  fitting.  Some  test  results  are 
presented  in  Section  5  and  conclusions  follow  in  Section  6. 


2  RATIONAL  APPROXIMATION  OF  RADI¬ 
ANCE  PROFILES 

Using  radiative  transfer  theory  one  can  relate  the  Planck  intensity  to  the 
upwelling  intensity  of  the  atmosphere.  This  lelationship  can  be  expressed  in 
terms  of  an  integral  equation  of  the  first  kind.  Specifically  the  relationship 
is  given  by 

Rip)=  r  B{p)wip/p)dp/p  (1) 

Jo 

where  yV{pfp)  is  a  kernel  weight  function  that  peaks  at  p  =  p  and  R{p) 
denotes  the  radiance  of  the  wavelength  channel  whose  weight  function  peaks 
o.t  p  =  p. 

If  we  set 

VV(z)  =  zc“*  and  s  =  i 

P 

then  equation  (1)  can  be  expressed  as  a  Laplace  transform  and  consequently 
we  can  solve  for  the  Planck  function 

We  can  think  of  B{p)  and  R{p)  as  transform  pairs.  Once  a  representation 
for  R  has  been  decided  upon,  then  B  can  be  determined  analytically  as  an 
inverse  Laplace  transform. 

_ The  /f-function  inversion  theory  of  Chandrasekhar,  (1950)  suggests  that 

R  should  be  represented  as  a  rational  function.  If  we  set 

R{p)  =  +  <^2P  H - +  dj+xp’  +  *•  .  (2) 

1  +  c,p 

then 

^(;)  _  5[l  .  <^2  dj+i  ^  wi 

s  s  ^  S2  ^^s  +  c< 

and 

J  t 

B(p)  =  di  d2P+ •••  + -^p^ +  (3) 

i=l 
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The  coefficients  in  equation  (2)  can  be  determined  from  the  data  by  ratio¬ 
nal  interpolation  or  by  nonlinear  least  squares.  Denote  the  radiance  value 
corresponding  to  p,  by  for  i  =  1,. . .  ,m  If  one  represents  the  radiance 
function  as  a  rational  function  of  the  form 


n,-x  aip’*  +  a2p"-* +  •••-»- a„+, 

°  -z, 


(n  -b  /  =  m  -  1) 


then  R  will  interpolate  (pi,R,)  if  and  only  if 


fllPi"  +  a2Pi"~*  + - !■  «n+l  +  Ri  d - b  =  Pi‘ Ri  (d) 


r 

The  coefficients  a,-  are  determined  by  solving  the  linear  System  (4).  The 
residue  form  of  R  can  be  computed  as  an  inverse  convolution  and  from  this 
one  can  determine  the  coefficients  for  equations  (2)  and  (3).  Generally  these 
equations  need  only  have  a  small  number  of  components.  It  is  possible  to 
get  a  good  fit  with  j  =  1,2  or  3  and  /  =  1  or  2.  If  one  takes  j  =  2  and 
/  =  2,  it  is  possible  that  for  some  data  sets  the  interpolating  function  will 
have  a  positive  pole,  (cj  <  0),  even  though  this  is  impossible  for  an  actual 
radiance  profile.  When  this  happens  the  weight  u),-  corresponding  to  the  pole 
wiU  be  much  smaller  than  the  weight  corresponding  to  the  other  rational 
component.  Thus  if  one  sets  Wi  =  0,  then  the  resulting  function  will  give  a 
good  approximation  to  the  data  points.  | 

In  practice  one  can  obtain  a  better  approximation  to  ^he  radiance  data 
and  avoid  the  problem  of  positive  poles  by  using  nonlinear  least  squares. 
The  rational  interpolating  function  can  be  taken  cis  startihg  approximation 
for  an  iteratively  computed  nonlinear  least  squares  fit  to  the  data  of  the 
form  ^ 

R{p)  =  xi X2P+ '■•  +  Xj+xp’ ^  (5) 

Initially  the  i,  coefficients  are  determined  from  the  coefficients  of  the  rational 
interpolating  function  of  the  form  (2).  Whenever  the  interpolating  function 
has  a  positive  pole,  i.e,  a  coefficient  c,  <  0,  one  sets  the  corresponding 
coefficient  ij+j+j+j  in  (5)  equal  to  1  and  the  weight  Xj+i+,-  can  be  set  equal 
to  0. 
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Figure  1.  Nonlinear  least  squares  radiance  function  for  set  of  TOYS  data 


The  choice  of  a  rational  function  to  represent  the  radiance  data  is  mo¬ 
tivated  by  the  physics  of  the  atmosphere.  In  practice  the  nonlinear  least 
squares  appro.ximation  gives  a  very  good  fit  to  the  radiance  data.  The  coef¬ 
ficients  of  the  rational  function  can  be  used  to  determine  a  Planck  function 
of  the  form  (3). 

In  Figure  I,  a  least  squares  fit  of  the  form  (-5)  is  given  for  a  typical  data 
set  obtained  from  the  NO.-VA  TIROS  Operational  Vertical  Sounder  (TOYS). 

3  GENERALIZED  EXPONENTIAL  INVERSION 

We  can  obtain  a  natural  generalization  of  our  previous  results  by  taking  a 
more  general  form  for  the  weight  function.  If  we  set 

W(2)  =  =  7*,re.xpi(2) 
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where 
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i'/‘r  (tfi) 

and 

expit{^)  =  exp(-2*/A:) 

then  it  follows  that 

R{p)=  r  B{p)w,{p/pA 

Jo  P 

If  R{p)  is  of  the  form  (2),  then  the  Planck  intensity  B{p)  =  Bk(p)  will 
be  of  the  form 

d  ' 

Bk{p)  =  di  +  djp  +  •  •  •  +  expjfc(c.p)  (6) 

i=i 

Note  that  equations  (3)  and  (6)  are  the  same  in  the  case  k  =  1. 

4  TOTAL  NONLINEAR  LEAST  SQUARES 

A  number  of  algorithms  have  been  developed  for  nonlinear  least  squares 
fitting.  Two  algorithms  using  gradient  methods  to  minimize  the  residual 
functions  were  developed  specifically  for  rational  fits  of  the  form  (5).  These 
algorithms  require  good  starting  approximations  and  may  converge  to  a  lo¬ 
cal  rather  than  a  global  minimum.  The  stamdard  Nelder-Meade  Simplex 
Algorithm,  however,  does  seem  to  give  satisiactory  results  even  though  con¬ 
vergence  is  rather  slow.  It  may  be  possible  to  speed  up  convergence  by 
taking  some  sort  combination  of  the  Nelder-Meade  and  gradient  methods. 

It  should  be  noted  that  the  standard  nonlinear  least  squares  techniques 
all  seek  to  minimize  the  residua!  errors  for  the  dependent  variable.  For  the 
type  of  radiance  fitting  problem  considered  here,  the  values  of  the  indepen¬ 
dent  variable  p  are  also  determined  from  satellite  data  and  may  involve  some 
error.  In  fact,  early  experiments  seem  to  indicate  that  the  computed  Planck 
function  is  more  sensitive  to  perturbations  in  the  pi  values  than  to  changes 
in  the  Ri  values.  In  order  to  reflect  this  in  the  data  fitting,  algorithms  have 
been  devised  for  total  nonlinear  least  squares  fits.  In  this  process  Ihe  resid¬ 
ual  is  defined  in  terms  of  the  sum  of  the  squares  of  the  distances  of  the  points 
(pi,  Ri)  to  tangent  lines  to  the  curve  (5)  at  the  points  (p,-,  i?(p,)).  Indeed  a 
parameter  t  has  been  incorporated  into  the  data  fitting  algorithms  to  specify 
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which  type  of  fit  is  used.  If  t  =  0,  then  the  usual  nonlinear  least  squares  fit 
is  used  to  minimize  the  residuals  of  the  dependent,  variable.  If  t  =  1,  then 
a  total  least  squares  fit  is  computed.  If  (  =  2,  then  a  rational  least  squares 
fit  is  computed  to  minimize  the  residuals  of  the  independent  variable.  For 
values  of  t  between  0  and  2,  intermediate  types  of  fits  are  computed. 


5  TEST  RESULTS 

Equation  (1)  is  an  integral  equation  of  the  first  kind  and  consequently  the 
problem  of  finding  a  solution  is  ill-posed.  In  addition  to  testing  the  optical 
measure  theory  method  we  also  tested  some  of  the  standard  methods  used 
for  solving  integral  integral  equations  of  the  first  kind. 

The  standard  techniques  involve  discretizing  the  equation  and  adding 
regularization  conditions.  The  problem  is  then  translated  into  a  linear  sys¬ 
tem  of  equations  with  quadratic  constraints.  The  constrained  problem  can 
be  solved  using  the  methods  proposed  by  Gander,  (1981).  However,  when 
these  techniques  were  applied  to  data  sets  obtained  from  the  NOAA  TIROS 
Operational  Vertical  Sounder  (TOYS),  the  resulting  matrices  had  very  low 
numerical  rank.  Consequently,  even  with  regularization  conditions  we  were 
unable  to  obtain  meaningful  results.  On  the  other  hand,  the  transform 
methods  have  smoothness  constraints  and  additional  structure  built  in.  The 
rational  form  (2)  can  be  computed  in  a  numerically  stable  manner  and  i,he 
coefficients  define  a  unique  Planck  function. 

The  interpolation  and  nonlinear  least  squares  algorithms  have  been  tested 
extensively  on  simulated  data  and  on  the  TOYS  data.  The  algorithms  are 
numerically  stable.  Figures  2  through  5  represent  Planck  functions  that 
were  generated  using  the  same  data  set  as  in  Figure  1.  Figure  2  shows  a 
semilog  plot  of  the  Planck  function.  The  pressures  are  given  in  millibars  on 
the  vertical  axis.  Figure  3  shows  the  Planck  function  for  the  data  set  de¬ 
rived  by  generalized  exponential  inversion  with  k  =  0.8.  Figure  4  shows  the 
Planck  function  derived  by  generalized  exponential  inversion  with  fc  =  1.2. 
Figure  5  shows  all  three  Planck  functions  plotted  as  functions  of  p  on  the 
same  axis  system. 
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The  question  that  still  is  outstanding  is  how  to  choose  j.  the  degree 
of  the  polynomial  component  of  the  ht,  and  how  to  choose  /,  the  i  umber 
of  hyperbolic  components.  Although  it  is  reasonable  that  both  of  these 
numbers  should  be  relatively  small,  it  is  not  possible  to  give  a  definitive 
answer  as  to  the  optimal  values  of  y  and  I  based  solely  on  the  TOYS  data  sets 
and  simulated  data.  It  is  anticipated  that  these  questions  will  be  answered 
when  future  data  sets  with  more  sensing  channels  in  the  appropriate  ranges 
become  available. 

6  CONCLUSIONS 

Although  transfer  theory  can  be  used  to  relate  the  Pianck  intensity  to  the 
upwelling  intensity  in  the  atmosphere,  the  relation  is  expressed  in  the  form 
of  an  integral  equation  of  the  first  kind.  Such  equations  are  ill-posed  and 
consequently  do  not  have  unique  numerical  solutions.  Adding  regularization 
conditions  does  not  solve  the  problem.  For  data  sets  like  those  obtained  from 
the  TOYS  data,  there  is  not  enough  information  to  determine  the  Planck 
function  using  matrix  methods  with  appropriate  smoothness  constraints. 
More  assumptions  relative  to  the  physics  of  the  atmosphere  must  be  added. 
In  this  regard  the  transform  methods  based  on  optical  measure  theory  seem 
to  work  well.  The  key  step  is  the  representation  of  the  radiance  profile  by  a 
rational  function.  This  can  be  accomplished  in  a  numerically  stable  manner 
using  nonlinear  least  squares  fitting  algorithmns. 
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